import numpy as np
import pandas as pd
import time
from datetime import datetime
from datetime import date
import scipy.stats as stats
from scipy.stats import norm
from scipy.optimize import minimize
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
from matplotlib import pyplot as plt

# 备注： 由于datetime.timedelta函数故障，故在计算计息天数时使用的月份*30方法可能造成一定误差，敬请谅解
# 假设： 假设名义本金为1个单位。


def st_simu(r, sigma, s_0, trade_day, kin_ob, seed=1):
    # 标的资产价格模拟路径
    s_t = [s_0]
    z = norm.rvs(0, 1, len(kin_ob), random_state=seed)  # 生成标准正态分布随机数
    s_temp = s_0
    for i in range(1, len(kin_ob)):
        temp = np.exp((r - sigma ** 2 / 2) * 1 / trade_day + sigma * np.sqrt(1 / trade_day) * z[i])
        s_temp *= temp
        s_t.append(s_temp)  # 生成对应日期的标的价格

    s_t = pd.DataFrame(s_t, index=kin_ob)
    return s_t


def state_judge(s_t, kin_ob, kout_ob, k_in, k_out):
    # 判断敲入/敲出
    out_flag = 0
    in_flag = 0
    out_time = []
    in_time = []
    for i in range(1, len(kout_ob)):
        if s_t.loc[kout_ob[i]][0] > k_out[i]:  # 判断敲出
            out_flag = 1
            out_time = kout_ob[i]
            break
    if out_flag == 1:
        return ["敲出", out_time]

    for i in range(1, len(kin_ob)):
        if s_t.loc[kin_ob[i]][0] < k_in:  # 判断敲入
            in_flag = 1
            in_time = datetime.strftime(kin_ob[i], '%Y-%m-%d')
            break
    if in_flag == 1:
        return ["敲入", in_time]
    else:
        return ["未敲出且未敲入", kout_ob[-1]]


def snow_price(r, sigma, trade_day, T, num_path, s_0, Y, NP, kin_ob, kout_ob, k_in, k_out):
    # 计算雪球期权价格
    snow_payoff = []  # 一份标的对应收益
    st_simu_df = pd.DataFrame()  # 储存模拟路径
    # 生成股票路径
    seed = np.linspace(1, num_path, num_path)
    state = []  # 状态
    knock_time = []  # 敲入/敲出时间
    knock_price = []
    for i in range(num_path):
        s_t = st_simu(r, sigma, s_0, trade_day, kin_ob, int(seed[i]))
        st_simu_df[str(seed[i])] = s_t.iloc[:, 0]
        res = state_judge(s_t, kin_ob, kout_ob, k_in, k_out)
        state.append(res[0])
        knock_time.append(res[1])
        knock_price.append(s_t.loc[res[1]][0])
        if res[0] == "敲出":
            out_time_list = res[1].split("-")
            k_count = (int(out_time_list[0])-2021) * 365 + (int(out_time_list[1])-1) * 30 + (int(out_time_list[2])-5)  # 计息天数
            snow_payoff.append(Y * NP * k_count / 365)
        elif res[0] == "敲入":
            snow_payoff.append(min(s_t.iloc[-1, 0] / s_0 - 1, 0))
        else:
            snow_payoff.append(Y * NP * T / 365)

    pricing_process = pd.DataFrame()  # 储存敲入敲出判断结果
    pricing_process['seed'] = seed
    pricing_process['state'] = state
    pricing_process['knock_time'] = knock_time
    pricing_process['knock_price'] = knock_price
    pricing_process['snow_payoff'] = snow_payoff
    print(pricing_process)

    t = (kin_ob[-1] - kin_ob[0]).days  # 剩余贴现时间长度T-t
    price = np.mean([snow_payoff[i] / ((1+r) ** (t/365)) for i in range(len(snow_payoff))])
    return price


def main():
    NP = 1  # 名义本金（自定义）
    s_0 = NP * 0.095  # 标的初始价格（融券成本, NP * q）
    sigma = 0.14  # 波动率
    r = 0.03  # 无风险利率
    Y = 0.16  # 敲出票息（年）/红利票息（年）
    T = 24 / 12 * 365  # 期限 （日）
    start_date = datetime.strptime('2021-01-05', '%Y-%m-%d').date()  # 开始日期
    end_date = datetime.strptime('2023-01-05', '%Y-%m-%d').date()  # 结束日期
    kin_ob = pd.bdate_range(start_date, end_date)  # 敲入观察日（工作日 暂时无法获取交易日）
    kout_ob = ['2021-04-06', '2021-05-06', '2021-06-07', '2021-07-05', '2021-08-05', '2021-09-06', '2021-10-08',
               '2021-11-05', '2021-12-06', '2022-01-05', '2022-02-07', '2022-03-07', '2022-04-06', '2022-05-05',
               '2022-06-06', '2022-07-05', '2022-08-05', '2022-09-05', '2022-10-10', '2022-11-07', '2022-12-05',
               '2023-01-05']  # 敲出观察日
    k_in = s_0 * 0.75
    k_out = [1, 1, 1, 0.995, 0.99, 0.985, 0.98, 0.975, 0.97, 0.965, 0.96, 0.955, 0.95, 0.945,
             0.94, 0.935, 0.93, 0.925, 0.92, 0.915, 0.91, 0.905]
    k_out = [element * s_0 for element in k_out]  # 敲出价格
    trade_day = len(kin_ob)
    num_path = 10000  # 模拟生成的样本路径数量
    prices = snow_price(r, sigma, trade_day, T, num_path, s_0, Y, NP, kin_ob, kout_ob, k_in, k_out)
    print('雪球期权价格约为本金的 %2.5f 倍。' % prices)


main()
